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I .  INTRODUCTION 


The  application  of  digital  filter  techniques  has  had  a  substantial 
impact  on  numerical  analysis  of  test  data  in  recent  years.  In  a  previous 
report,  reference  (1),  the  author  described  software  available  on  the 
BRL  CYBER  system  for  the  design  of  multiple  band  digital  filters  and 
differentiators.  It  is  the  purpose  of  this  report  to  illustrate  the 
techniques  by  which  such  filters  may  be  exploited  in  the  analysis  of  time 
series  data  taken  during  various  types  of  ballistic  experiments.  This 
report  will  also  describe  numerical  methods  for  verification  of  results 
obtained  through  the  use  of  digital  filters. 

It  is  important  to  define  carefully  the  type  of  data  to  which  these 
techniques  apply.  In  a  general  sense,  any  analog  data  taken  from  a 
system  which  has  governing  differential  equations  with  solutions  in  terms 
of  Fourier  series  is  a  candidate  for  such  analysis.  More  specifically, 
these  techniques  assume  the  data  to  consist  of  a  discrete  set  of  points 
sampled  from  an  analog  signal  at  equal  time  intervals.  In  addition,  it 
will  be  assumed  that  sufficient  care  has  been  taken  in  the  sampling 
process  to  avoid  aliasing,  (i.e.,  the  data  has  not  been  biased  by  over 
or  under-sampling)  and  that  no  component  of  the  recording  or  playback 
system  has  been  ’’overdriven"  in  frequency  response  or  amplitude. 

Finally,  it  should  be  noted  that  there  will  be  instances  in  which 
all  of  the  above  assumptions  are  satisfied,  and  vet  Fourier  analysis 
fails  to  provide  the  desired  information.  This  report  contains  at  least 
one  such  example.  The  important  point  to  be  remembered  is  that  these 
techniques  cannot  be  used  as  "black  box"  processes.  Moreover,  it  is 
possible  to  obtain  seemingly  reasonable  results  which  are  in  fact  totally 
erroneous  if  insufficient  care  is  given  to  the  analysis  process. 

II.  FREQUENCY  ANALYSIS  VIA  THE  FAST  FOURIER  TRANSFORM 

It  is  not  the  purpose  of  this  report  to  relate  in  detail  the  theory 
of  the  Fast  Fourier  Transform  (FFT) .  The  reader  will  be  assumed  to  have 
some  knowledge  of  the  subject,  and  is  referred  to  various  texts,  references 
(2,  5),  for  additional  information.  It  suffices  to  say  that  FFT  subroutines 
are  readily  available,  and  provide  the  background  for  this  report,  since 


J.IL  Walbert,  "Computer  Algorithms  for  the  Design  and  Implementation 
of  Linear  Phase  Finite  Impulse  Response  Digital  Filters",  BRL  lechnical 
Report  ( to  be  published) . 

p 

Bede-Liu,  editor,  Biyyital  Filter^  and  the  Fast  Fourier  Transforn, 
Halsted  Press,  1975. 

3 

L.R.  P.abiner,  B.  Gold,  Theory  and  Application  o£  Digital  Signal 
Prooessin£j  Prentice-Hall,  1975. 
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Fourier  analysis  without  the  FFT  is  far  too  time-consuming  computation¬ 
ally  to  be  feasible.  One  additional  comment  should  be  made:  in 
analyzing  ballistic  data,  one  deals  generally  with  relatively  small 
numbers  of  datij  points  (say  1000  to  2000)  during  a  specific  event.  For 
this  reason,  an  FFT  routine  based  on  powers  of  2  is  not  as  useful  as 
one  which  can  transform  a  data  set  of  arbitrary  length.  The  subroutine 
used  for  the  analysis  presented  in  this  report  is  one  given  to  the 
author  by  Dr.  Carl  de  Boor  at  the  Mathematics  Research  Center,  Madison, 
Wisconsin*.  A  listing  of  this  subroutine  in  BASIC  is  given  in  Appendix 
A. 

There  are  numerous  methods  for  using  the  FFT  to  compute  a  frequency 
spectrum  of  a  time  series,  and  perhaps  even  a  greater  number  of  technique 
by  which  one  may  emphasize  certain  spectral  characteristics.  References 
4,5,  and  b  provide  a  great  deal  of  insight  into  the  various  methods  which 
can  be  used.  For  the  processes  used  in  this  analysis  of  the  data 
discussed  in  this  report,  it  is  not  necessary  to  obtain  the  most  precise 
spectral  estimates  possible,  either  in  frequency  or  amplitude.  One  need 
only  be  able  to  determine  the  frequency  ranges  within  which  certain 
phenomena  of  interest  occur  in  order  to  design  filters  to  isolate  them. 
Once  these  individual  components  are  separated,  a  more  precise  frequency 
analysis  on  each  component  of  the  original  data  can  be  performed,  if 
desired.  It  suffices,  therefore,  to  use  a  simple  amplitude  spectrum, 
normalized  either  to  the  spectral  peak  or  by  some  other  suitable  means. 

For  example,  if  x(t)  are  the  data,  sampled  at  intervals  equally 
spaced  in  time  t,  then  the  discrete  Fourier  transform  of  x(t),  X(f),  is 
a  complex  valued  function 

X (f ) =R (f ) +i I (f ) , 


where  i  =-l  and  R  and  I  are  real  valued  functions  of  the  discrete 
frequencies  f.  The  magnitude  of  X(f),  denoted  by  |X(f)|,  is  given  by 


^Private  communication  during  a  visit  by  the  author  in  May ,  1979. 

sj  f 

R.B.  BLaakman ,  J.W.  Tukey,  The  Measurement  o/  Power  Spectra ,  Dover1 

Publications ,  1958. 


r>K.  Yuen,  D.  Fraser,  Digital  Spectral  Analpsisj  Pitman  Publishing, 

1979. 


6 


J.S.  Bendat,  A.G.  Piersol,  Random  Datav  Anabasis  and  Measurement 
Procedures,  Wiley  Interscience,  1971.  — " 
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Let 


M=max| X(f) | , 

taken  over  all  f.  Then  at  each  frequency  f,  the  value  S(.f)  of  a 
normalized  amplitude  spectrum  is  defined  by 


S (f)  =  | X (f) | /M. 


This  admittedly  crude  spectrum  will  be  shown  to  be  entirely  adequate 
for  the  first  stage  of  data  analysis. 

The  remainder  of  this  report  will  demonstrate  the  process  of  spectral 
analysis  on  5  types  of  data  curves.  These  same  five  examples  will  be  used 
to  illustrate  the  application  of  digital  filters.  In  each  example,  the 
data  has  not  been  converted  to  engineering  units;  the  vertical  axis  has 
been  labeled  signal  amplitude  for  reference.  This  was  deliberately  done 
to  illustrate  the  fact  that  the  vertical  scaling  of  the  data  is  not 
important  in  the  process  of  the  analysis.  Moreover,  each  example  is 
taken  from  test  data;  since  the  analysis  of  this  data  will  appear  in 
other  reports,  a  discussion  here  concerning  magnitudes  of  signal  compon¬ 
ents  in  engineering  units  would  serve  no  purpose. 

III.  EXAMPLE  1:  THE  HALF  SINE  PULSE  ANALYSIS 

The  data  curve  in  this  example,  shown  in  Figure  1,  is  the  record 
of  the  response  of  a  force  gage  mounted  on  a  small  mass  as  this  mass 
impacts  an  essentially  infinite  mass.  One  would  expect  the  force  to  be 
a  half  sine  pulse,  as  indeed  it  is,  but  there  is  an  additional  sinusoidal 
component  present  in  the  signal.  It  is  the  purpose  of  the  present  analysis 
to  separate  the  two  components  in  the  data. 

The  first  step  is  the  computation  of  the  Fourier  spectrum,  shown  in 
Figure  2.  Some  preliminary  comments  concerning  the  appearance  of  a 
spectrum  are  in  order.  In  general  terms,  the  more  one-sided  the  time 
series,  (i.e.,  the  larger  the  absolute  value  of  its  mean)  the  larger  the 
magnitude  of  the  spectrum  near  zero  on  the  frequency  axis.  In  this 
example,  one  therefore  expects  to  see  the  large  zero  Hz  component.  In 
some  cases,  this  component  may  in  fact  be  so  large  that  it  obscures  other 
components  in  the  spectrum.  One  solution  frequently  mentioned  in  the 
literature  is  to  adjust  the  data  so  that  it  has  zero  mean  prior  to  comput¬ 
ing  the  spectrum.  This  technique  works  quite  well  when  the  amplitude  of 
the  data  is  nearly  symmetric  about  its  mean;  in  other  instances,  such  as 
the  present  example,  mean  adjustment  has  little  or  no  effect  on  the 
spectrum.  At  any  rate,  it  is  not  a  goal  of  this  analysis  to  determine 
the  precise  amplitude  spectrum;  to  identify  the  frequency  ranges  of  the 
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Figure  1.  The  Half  Sine  Pulse:  Raw  Data 


13 


Figure  2.  The  Half  Sine  Pulse:  Normalized  Amplitude  Spectrum 


signal  components  the  spectrum  of  Figure  2  is  entirely  adequate.  Moreover, 
for  analysis  of  trends  or  other  baseline  variations  in  ballistic  data, 
such  as  in  Example  4,  adjustment  of  the  data  for  zero  mean  may  actually 
obscure  the  phenomenon  which  is  to  be  detected. 

Returning  to  the  spectrum  of  Figure  2,  one  notes  the  peak  at  zero 
Hz,  and  the  side  lobes  steadily  decreasing  in  amplitude  out  to  about 
4.1  kHz.  This  portion  of  the  spectrum  is  the  classic  shape  generated  by 
a  pure  half  sine  pulse.  The  rounded  peak  centered  at  6  kHz  is  the 
component  of  interest  here.  It  represents  that  portion  of  the  data  which 
is  extraneous  to  the  half  sine  pulse.  One  feature  of  this  peak  which 
should  be  noted  here  is  that  it  is  not  sharp.  This  is  an  indication 
that  more  than  one  isolated  frequency  is  present  in  this  data  component. 

There  are  two  possible  causes  for  this  resolution  problem: 

1)  there  are  numerous  data  components  present  with  frequency 
responses  too  close  together  to  be  distinguishable  in  the 
spectrum. 

2)  there  is  one  data  component,  the  frequency  of  which  varies  over 
the  duration  of  the  spectrum. 

In  the  first  case,  the  data  was  not  sampled  at  a  sufficiently  high  rate; 
in  the  second  case,  the  spectral  duration  must  be  shortened  for  more 
precise  determination  of  frequency  content. 

From  Figure  1,  it  can  be  seen  that  there  is  probably  only  one 
component  which  must  therefore  have  a  time-varying  frequency.  At  this 
stage  of  the  analysis,  there  is  no  need  to  recompute  the  spectrum; 
application  of  a  low  pass  filter  with  a  cutoff  frequency  somewhere  in  the 
range  of  2  to  4  kHz  will  effectively  remove  the  higher  frequency  component. 
In  this  instance,  choosing  a  cutoff  frequency  of  2  kHz  with  a  transition 
band*  from  2  to  4  kHz,  the  resulting  filtered  curve  is  shown  in  Figure 
3.  Depending  on  the  reason  for  conducting  the  analysis,  this  curve  may 
be  entirely  satisfactory.  There  is,  however,  one  obvious  difference 
between  the  curve  in  Figure  3  and  that  anticipated  from  Figure  1:  in 
the  original  data  curve,  the  start  of  the  pulse  is  defined  by  a  rather 
sharp  corner;  whereas  in  the  filtered  curve,  this  corner  appears  rounded. 

If  this  portion  of  the  curve  is  important  to  the  analysis,  then  this 
distortion  is  not  acceptable. 

The  "rounding"  is  an  undesirable  side  effect  of  the  filtering  process; 
a  technique  for  eliminating  this  problem  is  the  main  purpose  of  the 
present  example.  To  begin  with,  one  must  understand  how  and  why  the 
problem  occurs.  It  is  vital  to  realize  that  only  digital  data  can  have 
sharp  corners;  no  analog  system  has  components  with  the  infinite  response 

*All  filters  used  in  this  report  were  designed  and  applied  using  the 
software  described  in  Reference  1.  The  reader  is  assumed  to  be  familiar 
with  this  material  or  with  filter  design  in  general. 
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Figure  3.  The  Half  Sine  Pulse  Low  Pass  Filtered  with  Cutoff  Frequency  2  KHz 


necessary  to  produce  such  output.  The  analog  systems  are  nevertheless 
presumed  to  have  a  sufficiently  rapid  response  to  record  the  event,  and 
any  useful  digital  filtering  technique  should  not  degrade  the  quality 
of  that  response.  The  implication  is,  then,  that  the  response  of  the 
digital  filter  is  at  fault. 

Those  familiar  with  this  type  of  problem  in  numerical  analysis  will 
realize  that  there  are  two  solutions:  use  more  data  points  in  the 
vicinity  of  the  "corner",  or  use  a  shorter  filter,  i.e.:  one  with  fewer 
coefficients.  Each  of  these  solutions  has  drawbacks.  The  first  solution 
requires  sampling  the  data  at  a  higher  rate,  which  may  be  physically  or 
economically  impossible  for  some  types  of  data.  One  might,  of  course, 
sample  at  a  higher  rate  only  in  the  vicinity  of  the  corner,  but  non¬ 
constant  sampling  rates  pose  a  myraid  of  analysis  problems  and  are  best 
avoided  altogether.  In  addition,  increasing  the  sampling  rate  increases 
the  bandwidth  of  the  data.  This  means  that,  as  a  percent  of  bandwidth, 
any  low  pass  filter  of  a  specified  cutoff  frequency  has  smaller  pass  band 
and  transition  band,  which  may  in  general  lead  to  a  degradation  of  the 
filter  amplitude  response.  The  second  solution,  shortening  the  filter, 
may  also  cause  degradation  of  amplitude  response.  In  this  case,  one  is 
required  to  design  a  filter  response  with  what  may  prove  to  be  an  inadequate 
number  of  points. 

There  is  a  third  option.  The  major  underlying  assumption  in  Fourier 
analysis  is  that  the  data  to  be  analyzed  is  periodic.  This  assumption 
is  violated  by  the  half  sine  component  of  the  raw  data  curve.  Since  the 
corner  point  in  the  raw  data  set  can  be  identified  easily  by  numerical 
search,  one  can  reflect  the  data  set  through  this  point,  to  obtain  a  curve 
which  is  more  periodic.  Such  a  reflection  is  shown  in  Figure  4.  Assuming 
no  data  of  interest  occurred  prior  to  the  reflection  point  no  information 
is  lost.  Using  the  same  low  pass  filter  as  before  on  this  reflected  curve, 
and  then  zeroing  out  the  reflected  portion  and  truncating  the  data  set 
to  its  original  length,  one  obtains  the  filtered  curve  shown  in  Figure  5. 

The  corner  is  now  preserved  to  the  same  degree  of  accuracy  with  which 
it  appeared  in  the  raw  data.  This  reflection  process,  or  odd  periodic 
extension,  is  trivial  to  implement  in  a  computer  and  avoids  the  problems 
encountered  with  changes  in  filter  design  or  data  sampling.  One  could 
treat  the  trailing  end  of  the  half  sine  pulse  in  a  similar  fashion 
although  it  will  subsequently  be  shown  that  this  may  be  inappropriate. 

Although  Figure  5  shows  that  the  extraneous  component  has  been 
completely  removed,  this  example  is  obviously  simple.  How  one  can  determine 
the  adequacy  of  the  process  in  a  more  complex  data  set  will  be  demonstrated 
in  other  examples,  but  one  such  method  can  be  shown  in  this  example.  The 
Fourier  spectrum  of  the  curve  in  Figure  5  is  shown  in  Figure  6.  In 
comparison  with  the  spectrum  in  Figure  2,  it  will  be  seen  that  the  spectra 
are  virtually  identical  out  to  about  4.1  kHz,  beyond  which  the  spectrum 
of  the  filtered  data  set  shows  no  significant  content,  as  desired.  Thus, 
post-filter  spectral  analysis  (with  the  same  normalization  as  the  original 
spectrum)  will  verify  not  only  that  the  appropriate  components  have  been 
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Figure  4.  The  Half  Sine  Pulse  Reflected 
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Figure  5.  The  Half  Sine  Pulse  Low  Pass  Filtered  Using  the  Reflection  Technique 
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Figure  6.  Normalized  Amplitude  Spectrum  of  the  Filtered  Half  Sine  Pulse 


removed,  but  also  to  what  extent,  if  any,  the  remaining  components  have 
been  altered. 

It  is  an  indisputable  truth  that  any  analysis  of  this  type  of  data 
must  be  done  with  an  understanding  of  the  physical  and  mechanical  aspects 
of  the  experiment.  In  the  present  example,  a  small  cylindrical  deform¬ 
able  mass  impacted  an  infinite  nondeformable  mass.  A  careful  examination 
of  Figures  1  and  5  with  the  aid  of  a  straight  edge  will  reveal  the  fact 
that  the  force  curve  did  not  return  to  zero  after  the  event.  This  is 
evidence  of  the  fact  (physically  measured)  that  the  small  mass  was 
deformed  during  the  impact.  If  one  does  a  numerical  point-by-point 
subtraction  of  the  curve  in  Figure  5  from  the  curve  in  Figure  1,  one 
obtains  the  curve  in  Figure  7.  This  component  is  obviously  the  sinusoidal 
shockwave,  decaying  in  both  amplitude  and  frequency,  generated  internally 
in  the  small  mass  by  the  impact. 

Two  parts  of  this  curve,  at  the  beginning  and  end  of  the  half  sine 
pulse,  seem  to  deviate  from  its  otherwise  clean  appearance.  This  curve 
is  quite  useful  in  analysis  of  the  material  properties  of  the  small  mass. 
The  envelope  of  the  curve  is  a  measure  of  the  damping;  the  frequency  decay, 
used  in  conjunction  with  the  physical  dimensions  of  the  mass,  is  a  measure 
of  the  velocity  of  the  shock  through  the  mass,  as  well  as  an  indication 
of  the  compression  during  impact.  The  deviation  at  the  beginning  of  the 
curve  may  have  been  caused  by  the  impacting  surfaces  not  being  parallel 
initially,  setting  off  more  than  one  shock  wave,  with  phase  differences. 

The  deviation  near  the  end  of  the  half  sine  pulse  may  be  due  to  the 
rounded  corner  in  Figure  5,  but  a  numerical  simulation  of  this  impact  by 
a  Duhamel  integral  produced  a  curve  with  an  initial  sharp  corner  and  a 
terminal  rounded  corner.  For  this  reason,  it  may  not  be  appropriate 
(or  possible)  to  obtain  a  sharp  terminal  corner. 

IV.  EXAMPLE  2:  DAMPED  VIBRATION  ANALYSIS 

The  data  curve  of  this  example,  shown  in  Figure  8,  is  from  a  force 
gage  mounted  on  the  elevation  link  of  a  weapon  system.  The  impulse  of 
recoil  sets  off  vibrations  in  the  link  which  damp  out  in  time.  It  is  of 
interest  to  isolate  vibration  modes  and,  if  possible,  the  actual  impulse. 
The  spectrum  of  the  data,  Figure  9,  shows  peaks  of  interest  at  about  9, 

24,  and  47  Hz.  There  is  also  a  poorly  defined  peak  at  5  or  6  Hz,  which 
suggests  the  impulse  component.  The  spectrum  near  9  Hz  shows  two  peaks, 
but  these  are  so  close  it  seems  doubtful  they  can  be  separated  using 
filters.  In  this  analysis,  they  will  be  treated  as  a  single  component. 

Using  a  band  pass  filter*  passing  frequencies  from  3  to  7  Hz,  one 
obtains  the  curve  shown  in  Figure  10.  This  is  the  impulse,  It  should 


* Throughout  this  report >  band  pass  filtering  is  accomplished  in  two 
stages 3  using  a  low  pass  and  a  high  pass  filter .  The  process  is 
described  in  detail  in  section  VIII. 
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Figure  7.  Shock  Wave  Component  Isolated  from  the  Half  Sine  Pulse 
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Figure  8.  Elevation  Link  Recoil  Force:  Raw  Data 
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Figure  9.  Recoil  Force:  Normalized  Amplitude  Spectrum 
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Figure  10.  Recoil  Force  Impulse 


be  noted  that  in  isolating  each  component  in  this  example,  the  reflection 
process  demonstrated  in  example  1  was  used.  This  provides  an  accurate 
common  reference  point  on  each  curve.  The  accuracy  of  this  particular 
data  curve  can  be  verified  by  the  fact  that  a  numerical  simulation  of 
the  weapon  system  recoil  predicted  an  impulse  on  the  elevation  link  which 
agrees  with  the  curve  to  within  1%.  Applying  a  band  pass  filter  passing 
frequencies  from  8  to  12  Hz  produces  the  output  shown  in  Figure  11. 

This  represents  the  first  mode  of  vibration  and  shows  clearly  the  mode 
shape  and  the  damping. 

The  output  of  a  band  pass  filter  passing  frequencies  from  22  to  25 
Hz  is  shown  in  Figure  12.  This  curve  represents  a  higher  vibration  mode, 
with  at  least  one  phase  change  early  in  the  record.  The  particular  item 
instrumented  in  this  case  has  numerous  mechanical  stops  to  prevent  free 
vibrations;  these  stops  produce  the  phase  shifts  in  the  higher  modes. 

A  band  pass  filter  passing  frequencies  from  44  to  50  Hz  has  the  output 
shown  in  Figure  13.  This  mode  shows  at  least  two  phase  changes.  One 
feature  of  the  modes  shown  in  Figures  11-13  is  the  symmetry  about  the 
zero  amplitude  level.  This  symmetry  is  one  indication  that  these  modes 
have  been  isolated  correctly  by  the  filtering  process. 

There  is  another  method  by  which  one  can  determine  whether  all 
significant  portions  of  the  data  have  been  identified,  and  to  what  extent 
the  amplitudes  of  the  vibration  modes  are  correct.  This  method  is  simply 
to  form  a  point-by-point  sum  of  all  the  components.  Adding  the  curves 
in  Figures  10-13  produces  the  curve  of  Figure  14.  It  will  be  noted  that 
this  curve  bears  a  reasonable  resemblance  to  the  raw  data  curve  of  Figure 
8,  particularly  from  0.5  seconds  on.  It  is  clear,  however,  that  some 
significant  amplitudes  are  absent  prior  to  that  point.  A  re-examination 
of  the  spectrum  in  Figure  9  reveals  the  presence  of  other  peaks  which 
are  evidently  of  sufficient  significance  to  merit  further  study.  In 
particular,  the  peak  at  18  Hz  has  somewhat  higher  amplitude  than  the 
peak  due  to  the  impulse. 

A  band  pass  filter  passing  frequencies  from  15  to  20  Hz  produces 
the  output  of  Figure  15.  This  plot  shows  a  rapidly  decaying  signal,  with 
amplitude  concentrated  in  precisely  the  area  in  which  the  reconstructed 
curve  of  Figure  14  was  deficient.  Adding  this  curve  point-by-point  to 
that  of  Figure  14  produces  the  curve  of  Figure  16.  With  the  exception  of 
some -higher  frequency  content,  this  reconstructed  data  agrees  quite  well 
with  the  raw  data  in  Figure  8,  providing  good  confidence  in  the  results. 
One  might  continue  the  analysis  by  isolating  the  spectral  peaks  at  35, 

40,  and  53  Hz,  and  so  on.  It  is  clear  at  this  point,  however,  that  the 
major  components  have  been  identified. 

V.  CLOSED  COMPARTMENT  BLAST  PRESSURE  ANALYSIS 

As  was  indicated  earlier  in  this  report,  the  techniques  of  Fourier 
analysis  and  digital  filtering  are  not  always  completely  successful.  This 
next  example  illustrates  the  point.  Figure  17  is  a  plot  of  pressure  vs. 
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Figure  11.  Elevation  Link:  10  Hz  Vibration  Mode 
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Figure  12.  Elevation  Link:  24  Hz  Vibration  Mode 
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Figure  13.  Elevation  Link:  47  Hz  Vibration  Mode 
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Figure  14.  Sum  of  Impulse  and  10,  24,  and  47  Hz  Vibration  Modes 
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Figure  15.  Elevation  Link  Recoil  Force  18  Hz  Component 
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Figure  16.  Data  Curve  Reconstructed  from  Impulse  and  10,  18,  24,  and  47  Hz  Components 
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Figure  17.  Closed  Compartment  Blast  Pressure:  Raw  Data 


time  recorded  interior  to  a  closed  chamber.  The  chamber,  which  held 
several  containers  of  propellant  of  the  type  used  for  the  ballistic 
launch  of  projectiles,  was  subjected  to  penetration  by  a  shaped  charge 
jet;  this  in  turn  ignited  the  propellant.  The  data  therefore  contains 
a  pressure  component  due  to  the  burning  of  the  propellant,  a  blast 
component  with  probably  several  reflections  from  the  chamber  walls,  and 
quite  possibly  some  pressure  waves. 

Figure  18  is  the  spectrum  of  the  raw  data,  showing  the  main  non¬ 
periodic  contribution  at  zero  Hz  and  two  secondary  components  at  about 
450  Hz  and  1  kHz.  In  an  attempt  to  isolate  the  nonperiodic  zero  Hz 
component,  a  low  pass  filter  with  cutoff  frequency  of  200  Hz  is  applied 
to  the  data.  The  problem  which  arises  is  that  the  transition  band  of 
the  filter  is  not  sufficiently  narrow  to  eliminate  all  of  the  higher 
frequency  components.  As  a  result,  the  filtered  curve  appears  as  in 
Figure  19,  which  shows  the  pressure  component  from  the  burning  propellant 
with  some  additional  residual  higher  frequency  components. 

There  are  several  options  available.  First,  one  could  design  a 
filter  with  a  more  narrow  transition  band.  One  possible  consequence  is 
greater  deviation  from  the  design  specification  in  the  filter  pass  band, 
as  described  in  reference  1;  this  could  result  in  undesirable  amplitude 
distortion.  A  second  option  is  presented  by  the  fact  that,  having  low 
pass  filtered  the  data,  the  data  can  be  processed  at  an  increment.  For 
example,  the  bandwidth  in  the  spectrum  of  Figure  18  is  10  kHz.  If  one 
uses  every  other  point  (increment  by  2)  in  the  data  of  Figure  19,  then 
the  bandwidth  of  the  incremented  data  is  5  kHz.  The  significance  of  this 
fact  is  that,  since  the  filters  are  designed  based  on  fractions  of  band¬ 
width  rather  than  on  actual  frequencies,  the  range  between  the  200  Hz 
cutoff  point  and  the  450  Hz  spectral  component  is  now  twice  as  great. 

At  increment  4,  the  range  is  4  times  larger,  and  so  on.  Thus,  using  the 
same  filter,  the  transition  band  is  effectively  narrowed  by  incrementing 
the  data. 

While  this  second  option  works  in  this  example,  there  is  one  drawback 
incrementing  the  data  means  a  loss  of  data  points.  While  there  is  no 
information  lost  when  incrementing  is  done  in  conjunction  with  low  pass 
filtering,  one  can  no  longer  perform  point-by-point  operations  on  the 
incremented  data  with  some  other  component  at  a  different  increment, 
without  resorting  to  interpolation.  Depending  on  the  goal  of  the  analysis 
this  option  may  be  perfectly  acceptable. 

A  third  option  will  be  pursued  here,  to  demonstrate  another  useful 
technique  in  digital  data  analysis.  The  previous  example  demonstrated 
that  components  A,  B,  C,  and  D  of  a  raw  data  curve  R  could  be  summed 
point-by-point,  with  the  sum  being  nearly  equal  to  the  curve  R.  That  is, 


A+B+C+D  -  R 
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Figure  18.  Spectrum  of  Blast  Pressure  Data 
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Figure  19.  Output  of  a  200  Hz  Low  Pass  Filter  on  the  Blast  Pressure  Data 


It  would  seem  reasonable  to  expect  that  this  expression  can  be  manipulated 
so  that,  for  example, 


A  *  R-B-C-D 


That  is,  to  isolate  component  A,  one  could  just  as  well  isolate  components 
B,  C,  and  D,  and  subtract  them  point-by-point  from  the  original  data  R, 
leaving  an  approximation  to  A. 

This  method  works  especially  well  in  cases  like  the  present  one, 
where  there  is  a  large  amplitude  component  near  zero  Hz  with  numerous 
low  amplitude  components  close  by.  In  particular,  the  low  amplitude  peaks 
can  be  separated  out  one  at  a  time  by  high  pass  filtering,  beginning  with 
the  highest  frequency.  Using  this  method,  one  applies  a  high  pass  filter 
with  cutoff  frequency  of  950  Hz  to  the  raw  data.  A  copy  of  the  output 
at  this  stage  is  stored  for  future  use,  and  then,  the  output  is  low  pass 
filtered  at  a  cutoff  frequency  of  1.2  kHz.  The  output  of  this  2-stage 
band  pass  filtering  process  is  shown  in  Figure  20. 

This  curve  contains  several  features  of  interest.  The  initial  part 
of  the  curve,  the  first  three  positive  and  two  negative  peaks,  show  a 
nicely  decaying  wave.  Just  beyond  the  third  positive  peak,  this  pattern 
is  disturbed  in  both  phase  and  amplitude  by  the  second  of  the  two  main 
peaks  seen  in  the  raw  data  in  Figure  17.  Beyond  this  point,  the  curve 
has  little  recognizable  regularity.  The  important  fact  to  note  is  that 
the  clean  symmetry  of  some  of  the  components  in  the  previous  example  is 
absent  here. 

Continuing  the  analysis,  by  subtracting  the  output  of  the  highpass 
filter  with  cutoff  frequency  of  950  Hz  from  the  raw  data,  the  output  of 
Figure  21  is  obtained.  Comparing  Figure  21  with  Figure  19,  one  sees  that 
the  data  in  Figure  19  does  not  contain  the  full  450  Hz  component.  In 
order  to  separate  the  two  remaining  parts  of  the  data,  use  is  made  of  the 
incrementing  technique  described  earlier.  It  turns  out  that  by  incre¬ 
menting  the  data  in  Figure  21  by  two  (i.e.,  using  every  other  point)  the 
two  components  can  be  separated  using  a  high  pass  filter  with  a  cutoff 
frequency  of  200  Hz.  Since  the  original  data,  with  bandwidth  10  kHz, 
has  had  all  content  above  950  Hz  removed,  one  could  take  increments  as 
high  as  10  (a  bandwidth  of  1  kHz)  without  aliasing  the  data.  Using 
increment  2  simply  means  that  for  a  frequency  of  say  500  Hz,  one  has  20 
points  per  cycle  instead  of  the  original  40.  The  implication  is  that 
linear  interpolation  after  filtering  should  be  entirely  adequate. 

Incrementing  the  data  in  Figure  21  by  two,  applying  a  high  pass 
filter  with  cutoff  frequency  200  Hz,  and  using  linear  interpolation  on 
the  output  to  obtain  the  same  number  of  points  as  the  original  data  set, 
results  in  the  curve  in  Figure  22.  This  curve  also  shows  a  well  defined 
envelope  of  waves  up  to  the  end  of  the  third  negative  peak,  at  which 
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Figure  20.  Blast  Pressure  Data:  1  KHz  Component 
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Figure  21.  Blast  Pressure  Data  With  All  Components  Above  950  Hz  Removed 
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Figure  22.  Blast  Pressure  Data:  450  Hz  Component 
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Figure  23.  Blast  Pressure  Data:  Pressure  due  to  Propellant  Burning 


point  there  is  a  phase  change  and  a  loss  of  symmetry.  This  curve  bears 
some  resemblance  to  the  first  part  of  the  curve  in  Figure  20,  with  one 
major  difference;  in  Figure  22,  the  waves  seem  to  increase  in  frequency, 
rather  than  decrease  as  one  might  expect.  As  the  temperature  increases 
inside  the  chamber,  the  wave  velocity  (and  hence  frequency)  will  also 
increase;  this  may  in  fact  explain  the  increase  seen  here. 

At  any  rate,  subtracting  the  data  in  Figure  22  from  that  in  Figure 
21  results  in  the  curve  shown  in  Figure  23.  There  are  two  features  of 
this  curve  which  lead  one  to  question  its  accuracy.  First,  on  the  down¬ 
ward  slope  there  appears  to  be  some  residual  variation  from  the  450  Hz 
component.  More  importantly,  the  pressure  decreases  to  a  level  below 
zero  amplitude  for  a  short  time,  and  then  returns  to  zero.  This  negative 
pressure  level  is  really  a  part  of  the  blast  wave,  an  idealization  of 
which  is  shown  in  Figure  24.  The  data  of  this  example  contain  parts  of 
two  such  curves,  one  at  each  of  the  main  peaks  in  Figure  17.  Thus,  the 
analysis  here  has  failed  to  isolate  the  data  in  Figure  23  properly,  and 
the  actual  blast  curve  shape  is  partially  in  the  data  of  Figures  20,  22, 
and  23. 

The  curve  in  Figure  25  is  the  sum  of  the  curves  in  Figure  20,  22,  and 
23,  and  shows  that  the  mathematical,  if  not  physical,  components  have  been 
isolated  properly.  Thus,  this  particular  example  demonstrates  that  not 
all  physical  phenomena  are  amenable  to  Fourier  analysis.  The  above 
analysis  was  not  totally  in  vain,  however.  First,  the  data  in  Figure  23 
is  a  reasonable  representation  of  the  desired  curve.  Conceivably,  one 
could  use  a  filter  with  a  slightly  lower  cutoff  frequency  to  remove  the 
residual  higher  frequency  content  on  the  down  slope  of  this  data  curve. 

The  peak  value  would  seem  to  be  reliable,  as  well  as  the  maximum  value 
of  the  integral  of  this  pressure  curve.  Finally,  were  it  possible  to 
isolate  the  blast  pressure  curve  similar  to  that  of  Figure  24,  the 
integral  of  this  curve  should  approach  a  constant  value  which  is  a  measure 
of  the  force  applied  to  the  chamber  by  the  blast.  If  one  integrates  the 
450  Hz  component,  the  result  is  the  curve  shown  in  Figure  26.  Although 
there  are  some  early  cyclic  variations,  the  integral  does  finally  attain 
a  fairly  constant  value.  This  value  may  in  fact  be  a  measure  of  the 
blast  impulse.  It  should  also  be  noted  that  the  fact  that  the  final 
constant  value  of  this  integral  is  nonzero  is  proof  that  the  450  Hz 
component  does  not  consist  solely  of  pressure  waves. 

VI.  EXAMPLE  4:  ACCELEROMETER  BASELINE  CORRECTION 

A  frequent  problem  occurring  in  acceleration  data  is  an  apparent 
change  in  baseline  (zero  level)  during  the  event.  A  typical  example  is 
shown  in  Figure  27.  In  this  case,  the  level  change  is  obvious;  in  some 
cases  it  may  not  be  so  apparent.  There  is  a  way  of  checking  for  baseline 
variations  which  is  simple  yet  effective:  the  physical  component  being 
studied  has  a  velocity  which  returns  to  zero  after  the  data  event.  That 
is,  the  integral  of  the  acceleration  must  return  to  zero.  In  Figure  28, 
which  is  the  integral  of  the  data  in  Figure  27,  one  sees  that  there  is 
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Figure  24.  An  Idealized  Blast  Pressure  Curve 
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Figure  25.  Blast  Pressure  Data:  Sum  of  the  1  KHz,  450  Hz,  and  Non-periodic  Components 
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Figure  26.  Integral  of  450  Hz  Component 


(fl 

r 


oprn^dwy  i«u8vs 


45 


Figure  27.  Acceleration  Data  with  Baseline  Shift 
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Figure  28.  Integral  of  Uncorrected  Acceleration 


indeed  a  problem. 

In  analyzing  this  particular  difficulty,  one  commonly  made  mistake 
is  to  assume  that  at  some  point  during  the  data  event  an  abrupt  level 
change  occurred,  giving  the  baseline  the  appearance  of  a  step  function. 

While  this  is  possible  (although  the  change  can  occur  no  faster  than  the 
fastest  response  of  any  part  of  the  transducer/recording/playback  system) 
it  is  more  likely  that  the  analog  components  responded  in  a  ’’smooth” 
fashion.  Specifically,  if  a  step  function  such  as  that  in  Figure  29  is 
applied  to  the  input  of  an  analog  device,  an  output  of  the  form  of  that 
in  Figure  30  is  likely  to  occur.  Alternately,  it  may  be  that  a  ’’spike” 
will  occur  at  the  top  of  the  input  step,  as  seen  in  Figure  31.  This 
overshoot  may  result  in  an  analog  response  such  as  that  seen  in  Figure 
32.  The  point  is  that  one  should  anticipate  that  the  actual  baseline 
variation  has  the  appearance  of  the  curves  in  Figures  30  or  32,  or  some 
combination  thereof. 

A  further  indication  of  a  problem  in  the  raw  data  is  given  by  the 
presence  of  a  large  peak  at  zero  Hz  in  its  spectrum,  shown  in  Figure  33. 

The  spectral  contribution  of  this  component  includes  side  lobes  similar 
in  appearance  to  the  spectrum  of  Figure  6.  To  remove  the  component, 
a  low  pass  filter  with  a  cutoff  frequency  of  80  Hz  was  used.  This  cutoff 
frequency  is  well  below  the  frequency  of  the  first  vibration  mode  of  the 
system.  The  output  of  this  filter  is  shown  in  Figure  34.  The  similarity 
of  this  curve  to  that  of  Figure  32  is  readily  apparent.  If  this  curve 
is  subtracted  point-by-point  from  that  in  Figure  27,  the  result  is  the 
corrected  acceleration  data  shown  in  Figure  35. 

It  indeed  appears  that  the  level  shift  has  been  removed  by  this 
process,  and  in  fact,  the  spectrum  of  the  data,  Figure  36,  shows  the  zero 
Hz  component  and  its  side  lobes  to  be  absent.  The  real  proof  of  the  matter, 
however,  is  the  integral  of  the  corrected  acceleration  data,  shown  in 
Figure  37.  Whereas  the  final  value  of  the  integral  in  Figure  28  was  about 
75  units,  the  final  value  of  this  integral  is  indeed  zero  (precisely, 
the  computed  value  was  0.3  units,  although  had  the  integration  been  stopped 
earlier,  say  at  22  milliseconds,  the  computed  value  would  have  been  0.01 
units).  In  any  event,  this  integral  provides  immense  confidence  in  the 
baseline  correction  process.  Moreover,  the  shape  of  this  velocity  curve 
is  exactly  as  expected. 

VII.  EXAMPLE  5:  PRESSURE  WAVE  ANALYSIS 

This  final  example  demonstrates  the  fact  that  it  is  not  always 
possible  to  detect  the  presence  of  signal  components  by  looking  in  the 
time  domain.  Figure  38  is  a  plot  of  the  numerical  difference  between 
the  aft  and  fore  chamber  pressures  in  a  particular  gun  tube.  One 
phenomenon  frequently  present  in  this  pressure  difference  is  waves 
reflecting  between  the  base  of  the  projectile  and  the  rear  face  of  the 
tube.  In  this  particular  example,  it  would  appear  that  one  or  both  of 
the  pressure  curves  may  have  failed  to  correctly  represent  the  physical 
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Figure  29.  A  Step  Function 
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Figure  30.  Analog  Response  to  a  Step  Function 
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Figure  31.  Step  Function  with  Overshoot 
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Figure  32.  Analog  Response  to  Step  Function  with  Overshoot 
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Figure  33.  Spectrum  of  Uncorrected  Acceleration  Data 
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Figure  34.  Baseline  of  Acceleration  Data 
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Figure  35.  Corrected  Acceleration  Data 
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Figure  36.  Spectrum  of  Corrected  Acceleration  Data 
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Figure  37.  Integral  of  Corrected  Acceleration 
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Figure  38.  Pressure  Difference,  Aft  Gage  Minus  Forward  Gage 


phenomenon,  since  the  difference  curve  does  not  return  to  the  zero  level. 

The  spectrum  of  this  data,  Figure  39,  shows  the  large  zero  Hz  peak, 
reflecting  the  problem  noted  in  the  data.  There  are  well  defined  spectral 
peaks:  200  Hz,  420  Hz,  and  650  Hz.  Based  on  the  appearance  of  the 
original  data,  the  200  Hz  is  the  gross  cyclic  variation  in  the  curve. 

There  is  a  higher  frequency  present  in  the  data  which  is  not  noticeable 
in  the  spectrum.  The  420  and  650  Hz  components  are  not  apparent  in  the 
original  data. 

Nevertheless,  if  one  applies  a  band  pass  filter  to  the  data,  passing 
frequencies  from  400  to  600  Hz,  the  result  is  the  curve  shown  in  Figure 
40.  This  relatively  pure  waveform,  decaying  in  both  frequency  and 
amplitude,  is  a  typical  example  of  pressure  waves  found  behind  the 
projectile  in-bore.  What  is  especially  interesting  here  is  that  although 
the  original  record  was  of  questionable  validity,  the  wave  structure 
remained  intact. 

VIII.  SPECIAL  TECHNIQUES  FOR  DIGITAL  FILTER  APPLICATION 

In  the  foregoing  examples,  frequent  reference  has  been  made  to  the 
use  of  band  pass  filtering.  As  indicated  in  Example  2,  the  band  pass 
filters  employed  were  actually  two  filters,  a  low  pass  and  a  high  pass, 
applied  in  succession.  The  reason  for  using  a  two-stage  process  is 
simply  that  low  pass  and  high  pass  filters  with  the  required  frequency 
response  were  readily  available.  The  assumption  here  is,  of  course, 
that  it  is  accuracy  and  not  time  which  is  of  the  essence,  since  it 
takes  twice  as  long  to  apply  the  two  filters. 

This  process  is  shown  graphically  in  the  frequency  domain  in  the 
next  two  illustrations.  In  Figure  41  is  a  spectrum  containing  three 
peaks  (solid  lines) .  Suppose  it  is  desired  to  isolate  the  middle  spectral 
peak.  One  begins  by  applying  a  low  pass  filter  (dashed  lines)  with  a 
transition  band  situated  so  as  to  exclude  all  of  the  spectral  peak  on 
the  far  right  and  to  include  all  of  the  spectral  peak  in  the  middle. 

The  second  step,  shown  in  Figure  42,  is  to  apply  a  high  pass  filter 
(dashed  lines)  with  a  transition  band  situated  so  as  to  exclude  all  of 
the  spectral  peak  on  the  far  left  and  to  include  all  of  the  sole  remain¬ 
ing  spectral  peak  to  the  right.  Note  that  the  original  peak  on  the  far 
right  is  gone,  having  been  removed  by  step  1.  Obviously,  the  same  thing 
could  have  been  accomplished  by  using  the  high  pass  filter  first, 
followed  by  the  low  pass  filter. 

In  Figures  41  and  42,  the  spectral  peaks  were  sufficiently  far  apart 
to  allow  the  filter  transition  bands  to  fit  between  them.  Suppose  this 
is  not  true.  That  is,  suppose  the  spectrum  of  a  particular  data  set  has 
peaks  too  close  to  allow  the  transition  band  of  a  given  filter  to  fit 
between  them.  In  particular,  the  spectrum  of  Figure  43  has  four  peaks 
(solid  lines),  and  it  is  desired  to  remove  the  two  peaks  at  the  right. 
Using  the  same  low  pass  filter  as  in  Figure  41  (dashed  line),  it  will  be 
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Figure  39.  Spectrum  of  Pressure  Difference 
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Figure  40.  Pressure  Waves  Embedded  in  Pressure  Difference 
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Figure  41.  Band  Pass  Filter:  Step  1 
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Figure  42.  Band  Pass  Filter:  Step  2 
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Figure  43.  Multiple  Applications  of  the  Same  Filter 


seen  that  the  peak  on  the  far  right  is  removed,  but  the  second  peak  from 
the  right  is  not  completely  removed.  It  will  instead  be  reduced  in 
magnitude  by  about  half. 

One  option  is  to  design  a  new  low  pass  filter  with  a  more  narrow 
transition  band.  Another  option  which  is  very  convenient  is  simply  to 
apply  the  same  low  pass  filter  a  second  time.  This  has  the  effect  of 
squaring  the  magnitude  response.  In  the  pass  bands,  where  the  magnitude 
is  nearly  1,  the  squared  response  will  remain  nearly  1  (the  ripple,  or 
deviation,  will  of  course  be  squared,  but  if  it  is  small  to  begin  with, 
this  is  of  no  consequence).  In  the  stop  bands,  similarly,  the  values 
are  nearly  0  and,  when  squared,  remain  nearly  0.  It  is  in  the  transition 
bands  where  the  significance  of  the  procedure  becomes  apparent.  Positive 
numbers  less  than  1,  when  squared,  decrease  in  value.  The  net  effect  of 
two  filter  applications,  or  squaring  the  magnitude  response,  is  that  the 
transition  band  is  more  narrow,  as  shown  by  the  broken  line  (not  drawn  to 
precise  scale)  in  Figure  43. 

Thus,  if  the  second  peak  from  the  right  is  multiplied  by  G.5  in..each 
application  of  the  filter,  the  net  effect  of  two  applications  is  to  leave 
the  peak  with  0.5x0.5=0.25  times  its  original  amplitude. 

In  general,  a  filter  may  be  applied  as  many  times  as  desired, 
provided  the  amplified  deviation  in  the  pass  bands  can  be  tolerated.  If 
the  spectrum  of  the  data  is  sufficiently  clean,  even  large  amplitude 
deviations  in  the  pass  band  of  a  filter  can  be  corrected.  For  example, 
suppose,  as  in  Figure  44,  a  relatively  narrow  spectral  peak  (solid  lines) 
has  been  isolated  by  repeated  application  of  a  low  pass  filter  (dashed 
lines).  The  resulting  excessive  deviation  is  shown  by  the  arrows.  The 
output  of  this  filter  process,  being  a  single  frequency,  can  then  be 
multiplied  in  the  time  domain  by  the  amount  necessary  to  offset  the 
deviation.  This  is,  in  effect,  a  digital  amplifier.  It  must  be  emphasized 
that  this  works  only  at  isolated  frequencies  and  not  across  a  band  of 
frequencies . 

The  precise  amplification  factor  is  relatively  simple  to  determine. 
Construct  a  time  series  representation  of  a  sine  wave  with  amplitude 
varying  between  +1  and  -1  at  the  same  sampling  rate  as  the  original  data 
and  with  the  same  frequency  as  the  isolated  component.  Apply  the  digital 
filter  to  this  sinusoidal  data;  the  maximum  positive  output  amplitude  is 
the  inverse  of  the  amplification  factor.  For  example,  if  the  output 
sinusoid  has  maximum  positive  amplitude  0.9,  then  the  data  component 
must  be  multiplied  by  1/.9  to  correct  for  its  amplitude  distortion. 

In  cases  such  as  in  Example  1,  where  the  frequency  of  a  spectral 
component  varies  during  the  event,  this  amplification  correction  is  of 
no  value  since  the  spectral  peak  is  likely  to  be  wide  with  respect  to 
the  period  of  deviation  in  the  filter  pass  band  as  shown  in  Figure  45. 

It  is  clear  that  no  single  amplification  factor  will  correct  the  distor¬ 
tion  resulting  here.  It  should  also  be  clear  that,  if  there  is  more  than 
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Figure  44.  Amplitude  Correction  for  Excessive  Pass  Band  Deviation 


Figure  45.  Noncorrectable  Amplitude  Distortion  (I) 


one  distinct  frequency  that  has  been  isolated,  the  distortion  may  not 
be  correctable.  In  Figure  46,  for  example,  each  peak  will  be  distorted 
by  a  different  amount;  no  single  correction  factor  can  be  applied  to 
the  output  to  correct  this  distortion. 

IX.  NUMERICAL  DIFFERENTIATION  BY  DIGITAL  FILTER 

Numerical  differentiation  has  always  been  at  best  a  tenuous  process. 
The  primary  difficulty  with  most  numerical  schemes  is  their  inherent 
instability  at  high  frequencies.  In  this  section  it  will  be  shown  that 
digital  filters  designed  as  differentiators  are  stable  and  highly 
accurate.  The  theoretical  basis  for  the  design  of  a  differentiating 
filter  is  quite  simple.  Let  f(t)  denote  a  differentiable  function  of 
time  t,  and  suppose  F(w)  is  the  Fourier  transform  of  f,  where  w  is  the 
circular  frequency. 

Then 


f(t) 


e1U)tF(m)dm 


Similarly,  if  G(io)  is  the  Fourier  transform  of  the  derivative  of  f(t) 
with  respect  to  t,  f'(t),  then 


f'(t) 


c» 

■  */ 


ela3tG(to)d{ja 


From  the  first  of  these  two  equations, 


f'(t) 


ei^t  (iojF  taO  )  da), 


so  that 

G  (co)  =  ia)F(a)). 


It  follows,  then,  that  differentiation  in  the  time  domain  corresponds 
in  the  frequency  domain  to  multiplication  by  w  and  positive  rotation  of 
the  phase  angle  through  tt/2  radians.  In  particular,  the  amplitude 
response  is  just  |w|,  since  |i|=l.  Now,  the  variable  uj  is  the  circular 
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Figure  46.  Noncorrectable  Amplitude  Distortion  (II) 


frequency,  expressed  in  radians  per  second.  One  could  write 


0)  =  27lf, 


where  f  is  the  frequency  in  units  of  cycles  per  second,  or  Hz. 

In  the  case  of  the  discrete  Fourier  transform,  the  Nth  discrete 
frequency  f^  is  given  by 


fN  =  (N-l) /  [M(At)  ]  , 


where  M  is  the  number  of  time  series  points  spaced  At  apart.  Then 


and  (N-l) /M  is  the  fraction  of  the  sampling  frequency.  The  amplitude 
response  of  a  differentiating  filter  can  then  be  represented  as  in 
Figure  47.  It  is  an  approximation  to  the  straight  line  through  the 
origin  and  through  the  point  {h,h) ,  where  the  first  coordinate  is 
amplitude  and  the  second  coordinate  is  fraction  of  sampling  frequency. 

After  applying  this  differentiating  filter  to  a  data  set,  the  differen¬ 
tiated  data  is  multiplied  by  2Tr/At  to  correct  the  amplitude. 

The  filter  response  of  Figure  47  is  a  good  approximation  to  the 
ideal  differentiator  out  to  about  80%  of  the  bandwidth.  Beyond  that 
point,  the  amplitude  falls  off  rapidly  to  0.  It  is  this  return  to  zero 
amplitude  that  provides  the  stability  at  high  frequencies.  Clearly,  the 
amplitudes  of  frequencies  above  80%  of  bandwidth  are  greatly  decreased. 
Instead  of  increasing  the  noise  content  of  the  differentiated  signal, 
the  high  frequencies  are  actually  suppressed.  Since  the  magnitude  response 
of  the  filter  is  constrained  by  design  to  be  odd  periodic,  it  must  return 
to  zero  at  each  end;  it  is  not  possible  to  design  a  stable  differentiator 
over  the  full  bandwidth. 

As  long  as  the  data  set  contains  no  useful  information  beyond  80% 
of  the  bandwidth,  this  differentiator  performs  admirably  well.  For 
example,  if  one  computes  by  the  trapezoidal  rule  the  integral  of  the 
data  in  Figure  17,  the  result  is  shown  in  Figure  48.  On  applying  the 
differentiating  filter  of  Figure  47  to  the  integral  in  Figure  48,  one 
obtains  the  curve  shown  in  Figure  49.  This  curve  is  virtually  identical 
to  that  of  Figure  17;  in  fact,  the  difference  between  the  two  curves  is 
shown  in  Figure  50.  Indeed,  the  spectrum  in  Figure  18  shows  that  there 
is  very  little  content  beyond  80%  of  the  bandwidth  (8kHz). 
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Figure  47.  Frequency  Response  of  a  Differentiating  Filter 
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Figure  48.  Integral  of  Blast  Pressure  Data 


Figure  49-  Derivative  of  Integral  of  Blast  Pressure  Data 
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Figure  50.  Difference:  Figure  17  Minus  Figure  49 


One  could  not  expect  the  differentiating  filter  to  do  as  well  if 
the  input  data  contained  significant  frequencies  beyond  the  80%  point. 

For  example,  the  acceleration  data  of  Example  4  was  deliberately 
sampled  sparsely  so  that  its  spectrum  (Figures  33  and  36)  has  significant 
content  out  to  the  very  end.  Therefore,  applying  the  differentiating 
filter  to  the  integral  in  Figure  37  results  in  the  curve  in  Figure  51. 
While  this  curve  is  similar  to  that  of  Figure  35,  it  is  obvious  that 
some  high  frequency  content  was  lost.  In  fact,  the  difference  curve 
of  Figure  52  has  nearly  the  same  maximum  and  minimum  amplitude  as  that 
of  Figure  51 . 

If  it  should  become  necessary,  in  the  course  of  analysis,  to 
differentiate  the  data  in  Figure  37,  there  is  a  straightforward  method 
by  which  this  may  be  accomplished.  If,  in  a  particular  time  series,  one 
interpolates  between  data  points  to  double  the  total  number  available, 
then  the  apparent  sampling  rate  is  doubled,  as  is  the  bandwidth  of  the 
data.  Using  a  simple  linear  interpolation  on  the  data  of  Figure  35,  its 
spectrum  now  is  as  shown  in  Figure  53.  This  spectrum  is  identical  to 
that  of  Figure  36  up  to  the  8  kHz  point.  The  ’’fabricated”  part  of  the 
data,  from  8  kHz  to  16  kHz,  shows  a  relative  symmetry  about  the  8  kHz 
point,  but  with  a  rapidly  decreasing  amplitude.  The  important  point  to 
note  is  the  absence  of  significant  spectral  content  beyond  80%  of  the 
new  bandwidth. 

Applying  the  differentiating  filter  to  the  integral  of  the  expanded 
acceleration  data,  and  then  incrementing  the  output  by  two  to  return  to 
the  original  bandwidth,  results  in  the  curve  shown  in  Figure  54.  This 
curve  is  much  closer  to  that  of  Figure  35  than  was  that  in  Figure  51. 

In  fact,  Figure  55  shows  the  difference  between  Figure  35  and  Figure  54 
to  be  insignificant  except  for  the  two  high  frequency  peaks,  one  positive 
and  one  negative,  at  the  center  of  the  actual  data  event.  While  this 
technique  has  been  relatively  successful,  it  does  demonstrate  that  using 
a  sampling  rate  with  only  two  samples  per  cycle  at  the  highest  frequency 
of  interest  may  be  enough  to  represent  the  data,  but  not  enough  to  analyze 
it. 


More  precisely,  sampling  rate  and  useful  frequency  content  can  be 
linked  as  follows:  at  the  100%  of  bandwidth  point  in  the  spectrum,  one 
has  two  samples  per  cycle,  while  at  the  50%  bandwidth  point,  one  has 
four  samples  per  cycle.  It  follows  that  at  80%  of  the  bandwidth,  one 
has  about  2.5  samples  per  cycle  (although  the  relationship  is  nonlinear). 
That  is,  a  minimum  reasonable  data  sampling  rate  is  3  samples  per  cycle 
at  the  highest  frequency  of  interest.  In  actual  practice,  especially 
for  analysis  purposes,  it  is  not  desirable  to  have  significant  spectral 
content  too  close  to  the  80%  point,  and  hence  a  much  more  useful 
minimum  sampling  rate  is  4  or  even  5  samples  per  cycle.  Moreoyer,  it 
has  been  shown  that  components,  the  frequency  of  which  changes  through 
the  data  event,  can  be  isolated  from  the  data  for  further  study.  If  this 
additional  analysis  includes  advanced  Fourier  techniques  such  as  over¬ 
lapping  transforms  for  precise  frequency  determination,  8  to  10.  samples 
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Figure  51.  Derivative  of  Integral  of  Acceleration  Data 
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Figure  52.  Difference:  Figure  51  Minus  Figure  35 


77 


Figure  53.  Spectrum  of  Expanded  Acceleration  Data 


I 


aprnnduiy  ^uBjs 


I 

i 

i 


78 


Figure  54.  Derivative  of  Integral  of  Expanded  Acceleration  Data 
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Figure  55.  Difference:  Figure  35  Minus  Figure  54 


per  cycle  at  these  frequencies  may  be  necessary. 

To  further  complicate  the  picture,  it  may  not  be  possible  to 
determine  an  appropriate  sampling  rate  prior  to  conducting  an  experiment, 
since  some  or  many  of  the  frequency  ranges  may  be  unknown.  In  this 
respect,  sampling  rates  should  be  tailored  primarily  to  the  frequency 
responses  of  the  various  record/playback  devices  and  to  the  capabilities 
of  the  transducers  involved.  Total  system  frequency  response  cannot  be 
increased  by  higher  sampling  rates,  but  data  can  certainly  be  lost  by 
sampling  rates  which  are  too  low.  In  short,  instrumentation  frequency 
responses  and  sampling  rates  must  be  chosen  in  consideration  of  the 
overall  analytical  plan  in  order  to  provide  truly  useful  data.  References 
7,  8,  and  9  contain  more  detailed  discussions  of  these  concepts. 

X.  CONCLUSIONS 

In  the  analysis  of  ballistic  data  such  as  pressure,  strain,  and 
acceleration,  the  Fourier  spectrum  is  an  extremely  useful  tool.  For 
nonreal-time  data  analysis,  various  numerical  manipulations  can  be 
performed  to  enhance  the  applicability  of  Fourier  techniques.  These 
manipulations  include  odd  periodic  continuation,  linear  interpolation, 
and  multiple-pass  digital  filter  applications.  While  several  less  than 
successful  examples  have  been  presented,  one  can  nonetheless  conclude 
that  Fourier  analysis  is  applicable  even  in  the  presence  of  nonstationary 
frequencies.  It  can  also  be  used  to  isolate  aperiodic  phenomena  such  as 
baseline  variations. 

Numerous  methods  have  been  presented  by  which  the  validity  of  the 
results  of  Fourier  analysis  can  be  verified,  and  by  which  amplitude 
errors  can,  in  some  cases,  be  corrected.  It  must  be  emphasized  that  in 
any  such  application,  a  great  deal  of  engineering  judgment  is  required 
for  successful  analysis.  There  is  no  substitute  for  thorough  knowledge 
of  engineering  and  physical  principles  as  well  as  mathematical  theory. 
Finally,  specific  results  concerning  the  mechanical  behavior  of  a  system 
or  revealing  some  fundamental  physical  or  chemical  property  should  be 
repeatable.  The  Fourier  spectrum  is  a  statistic;  specific  conclusions 
based  on  this  type  of  analysis  must  be  made  from  a  sufficiently  large 
sample  of  data. 


y 

R.K.  Otnes,  L.  Enochson ,  Digital  Time  Series.  knalgsisj  Wiley  Inter- 
science,  1972. 

g 

R.  W.  Hamming,  Numerical  Methods  for  Scientists  and  Engineers,  McGraw- 
Hill,  1962. 


9 

A.  Papoulis,  The  Fourier  Integral  and  its  Applications,  McGraw-Hill, 
1962. 
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APPENDIX  A 


A  COMPUTER  SUBROUTINE  FOR  THE  FFT 
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580  Rng  1  e=Twop  i  •••'Th  i  r d 
5  9  0  0  m  e  q  ar  =  C  0  S  <  R  n  g  1  e  > 
600  Omegai =-S I N ( Rng  1  e 


610  fir  g  r  - 1 
620  firg i =0 

630  Sec  ondsBef ore*Rf t er 

640  F i rst sNoumone*Sec ond 

650  FOR  Index3add=0  TO  Nowmone*flf ter  STEF'  After 
668  FOR  I a= 1  TO  After 
670  I ndex 1 =F i rst  + 1  a 

63 0  FOR  I ndex3= I ndex3add+ 1  a  TO  I ndex3add+ I a+Fourt h  STEP  Third 

6  9  0  V  a  1  u  e  r = 2  i  n  <  1 ,  I  n  d  e  x  1  > 

700  V a 1 u  e i =  Z i n  < 2 , I  n  d  e  x  1  > 

710  FOR  I ndex2= I ndex 1 -Sec ond  TO  I ndex 1 -First  STEP  -Second 

7  20  Real  =Va1  u  e  r  *  fi  r  g  r  -  V  a  1  u  e  i  *  A  r  g  i  +  2  i  n  (1,1  n  d  e  :>■:  2  > 

730  Val uei =Val uer*F}rgi +  Val uei *firgr+2i n<2, Index2> 

740  Val uer=Real 

750  NEXT  I ndex2 
760  Zout ( 1 p I ndexS > =Val uer 
770  Zout <2, I ndex3 ) =  Val ue i 
780  I n  d e x 1  =  I n  d  e x 1 +  fi  f  t e r 

798  NEXT  Index3 

80O  Real  sArgr*Qriiegar~Argi  *Dmegai 

810  flrgi -firgi *Omegar +firgr *0megai 

820  firgr=Real 

830  NEXT  la 

840  NEXT  I ridexSadd 

850  SUBEND 
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